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Abstract 

Finite-temperature micromagnetics simulations are employed to study the magnetization- 
switching dynamics driven by a field applied at an angle to the long axis of an iron nanopillar. A 
bi-modal distribution in the switching times is observed, and evidence for two competing modes 
of magnetization-switching dynamics is presented. For the conditions studied here, temperature 
T = 20 K and the reversal field 3160 Oe at an angle of 75° to the long axis, approximately 70% 
of the switches involve unstable decay (no free-energy barrier) and 30% involve metastable de- 
cay (a free-energy barrier is crossed). The latter are indistinguishable from switches which are 
constrained to start at a metastable free-energy minimum. Competition between unstable and 
metastable decay could greatly complicate applications involving magnetization switches near the 
coercive field. 

PACS numbers: 75.75.+a, 75.60.Jk, 75.40.Mg, 85.70.Ay 
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I. INTRODUCTION 



Nanoscale magnetic devices play important roles in many applications, including sensor 
technology and magnetic recording. Fabricating such devices requires innovative techniques 
that rely on knowledge of the sensor's structure and dynamics at the nanoscale. The role 
of a micromagnetic simulation in the design process, consequently, is to represent such 
systems with the appropriate resolution, i.e., preserving the actual physical dynamics while 
maintaining a reasonable simulation effort. This has historically been achieved through 
the use of semi-classical equations that govern the motion of individual spins on a lattice 
mapped from a physical magnetic system. In our work, the Landau-Lifshitz-Gilbert (LLG) 
equation- provides such dynamics for the spins and includes finite-temperature effects by 
incorporating a stochastic field that obeys a fluctuation-dissipation relation. 

Here we model the magnetization switching of a system with strong uniaxial shape 
anisotropy: an iron nanopillar. When subjected to a magnetic field obliquely aligned with 
respect to the pillar axis, the pillar exhibits a bimodal distribution of switching times. 2 This 
is a feature which may have important ramifications for the application of such nanopillars in 
real-world devices, which typically rely on a single, consistent decay mode. The traditional 
picture of a single free-energy barrier which must be crossed, often used to describe the decay 
of a metastable state, appears to be insufficient for this system. Therefore, we analyze the 
bimodal switching behavior of simulated nanopillars using information from temporal phase 
portraits and also apply absorbing Markov- chain techniques^ 1 ^^^ to transition matrices 
obtained from the simulations. 

The rest of this paper is organized as follows. We first present our computational model 
and our implementation of the LLG equation in Section [Til Section II III is divided into three 
parts, which collectively discuss results obtained from the simulations. Section UlIAl explores 
the phase portraits of the energy during switching, while Section IIIIBI and Section IHI CI 
provide information about the free energy of the system based on analysis of transition 
matrices and projective dynamics, respectively. Our conclusions are presented in Section HVl 
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II. MODEL AND NUMERICAL METHOD 



Our numerical model is motivated by nanopillars fabricated by von Molnar and collab- 
orators.- 1 ^ Using scanning-tunneling-microscopy-assisted chemical vapor deposition, they 
constructed elongated iron nanoparticles, each approximately 10 x 10 x 150 nm 3 . Their 
results indicated that the field-driven magnetization switching in these particles at field 
strengths near the coercive field is initiated by localized nucleation, followed by thermal 
activation over a free-energy barrier.— illil&Ii This is expected for systems whose dimen- 
sions are large enough to support a non-uniform magnetization,—^ in contrast to the single 
coherent rotation mode assumed by smaller particles via the Stoner-Wohlfarth model. 17 

In order to numerically investigate these nanopillars, we use a coarse-grained, cuboid 
computational lattice in which each cell represents the net magnetization of the correspond- 
ing volume in the physical system. To ensure that the magnetization density is uniform at 
length scales below the cell volume, the lattice spacing is chosen smaller than the exchange 
length of 2.6 nm, obtained from the material properties of bulk iron.— This criterion yields 
a regular lattice with the dimensions 6 x 6 x 90 (N = 3240 spins, Aa; = 1.6667 nm), which 
has a single classical Heisenberg spin at each site. The time evolution of this spin, m(r\), is 
controlled by the LLG equation,-^^ 
drh(fi) _ 7q 
dt 1 + a 

which updates every site on the lattice at each computational time step. At each site i, the 
local field H(fl) determines the direction of change of the magnetization m(fl) during the 
next integration step. The parameters 70 = 1.76 x 10 7 Hz/Oe, m s = 1700 emu/cm 3 , and a = 
0.1 represent the electronic gyromagnetic ratio, the saturation magnetization of bulk iron, 
and a phenomenological damping parameter, respectively. Values for these parameters are 
consistent with the material properties of bulk iron, and are discussed in previous work.-^^i 
The total local field, H(ri), is a sum of the individual fields, which include the dipole field 
H D , the applied field H z , the exchange field H E , and the stochastic thermal field H T . 

Only nearest-neighbor interactions are considered for the exchange field H E , while the 
long-range dipolar interaction, H B , couples all computational cells to each other. Conse- 
quently, the dipolar term constitutes the largest computational task during an integration 
step of the simulation, scaling as 0(N 2 ) for a brute force calculation. This term is reduced 
to a linear dependence on the number of cells, N, by using the Fast Multipole Method,—^ 



-> 01 -> 
H(fi) x H{fi 



(1) 



resulting in simulation rates of approximately 0.005 ns/cpu-hour on a 3 GHz Intel Pentium 
4 processor with 512 MB of memory using the C++ Psimag library. 23 

The temperature is set to 20 K and is included in the simulation through the stochastic 
field H T , whose components are Gaussian distributed with mean zero and variance deter- 
mined by the fluctuation-dissipation relation, 



{Hj{r i: t)H^{f^')) = ^Ls^it-t'l (2) 

where fee is Boltzmann's constant, V is the volume of an individual computational cell, T 
is the absolute temperature, 5^ and <5g 7 are Kronecker deltas over the lattice sites i,j and 
directions /3, 7, respectively, and 5(t — t') is a Dirac delta function of the time difference, 
t — t'. This equation implies that the magnitude of the thermal field scales linearly with 



the square root o 
discussed in Ref. 



the temperature. The details of the integration of the stochastic field are 
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The simulation begins with the pillar in an applied field of magnitude 3160 Oe and angle 
75° with respect to the long axis of the pillar (Figfj^a)). Once equilibriated, the spins are 
relaxed in the direction of the initial applied field (FigJTJb)). The value of the applied field 
is then changed sinusoidally to a direction anti-parallel to the initial field over a time period 
of t = 0.125 ns, H(t) = H(0) cos(7rf/0.125 ns), t G [-0.125 ns, 0]. The ^-component of H(0) 
is negative (Figfjjc)) This is done to ensure that excess energy is not artificially pumped 
into the system. Upon completion of the field reversal, the simulation time is set to zero 
and the field is kept static at H(0) for the remainder of the trial. 



III. NUMERICAL RESULTS 

Under the conditions described in Section [Til 42 magnetization-switching simulations were 
performed. The switching time t B is defined as the first-passage time to Xli=ijv m 2( r l) < 0, 
with m z the component of the magnetization along the long axis of the pillar, and time mea- 
sured from the completion of the initial field reversal. From these 42 trials, the cumulative 
distribution of the switching times, shown in Fig. 12(a), indicates at least two characteristic 
time scales and the existence of more than a single switching path. About 70% of the sim- 
ulations (29 runs) switched almost immediately (t s < 2 ns), while the remaining 30% (13 
runs) exhibited switching times at least an order of magnitude larger. Since direct compar- 
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Start Equilibrated Field reversal Final 

t=0 equilibrium 



FIG. 1: (Color online.) Schematic of pillar magnetization (small arrows, black online) and the 
applied field (large arrows, red online) at different times during a trial, (a) Initially, at t = —0.25 ns, 
the magnetization is aligned with the long axis of the pillar and the applied field has a value of 
3160 Oe, oriented at 75° with respect to the long axis, (b) Some time later, at t = —0.125 ns, 
the magnetization is relaxed and the applied field begins its sinusoidal reversal as described in the 
text, (c) The field reversal is complete at t = ns. (d) The final equilibrium of the simulation 
after magnetization switching has occurred. 



ison of the magnetization of individual runs did not reveal any differences in the switching 
mechanism, the trials were divided into two groups based only on the observed switching 
times: "fast" decay (t s < 2 ns) and "slow" decay (t s > 2 ns). Below we show that the 
slow-mode statistics are the result of a process which must traverse a free-energy landscape 
characterized by a metastable well that the system must escape to reach the lowest avail- 
able free-energy state. This is accomplished by the collective effect of many random thermal 
fluctuations which can eventually cause the system to surmount the free-energy saddle point 
that separates the metastable well from the global free-energy minimum. 

Thermally activated barrier-crossing rates from processes such as this are described by 
an Arrhenius-like expression and a switching-time probability density function (pdf) which 
assumes an exponential form. To find the lifetime, r, of each mode, the mean switching time, 
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FIG. 2: (Color online.) Cumulative switching-time distribution for all 42 full runs (a) and a com- 
parison of quenched-relaxed and slow runs (b). It is clear from (a) that at least two characteristic 
time scales are present in the switching statistics for this system. From the measured lifetimes of 
the trials we obtain (t) = 1.2 ns with a standard deviation at = 0.4 ns for the fast mode (29 trials), 
and for the slow mode (t) = 20.9 ns, to = 7.8 ns, and r = at = 13.1 ns (13 trials). Also shown are 
the lifetimes for the QR trials with (t) = 14.2 ns, to = 1-6 ns, r = 12.6 ns and at = 10.8 ns. The 
fitted expressions, D(t) = rj{t — to)(l — exp(— (t — to)/r)), where r\ is the Heaviside step function, 
are also shown in (b). 

(t) = (1/iV) Yli^ii an d the standard deviation, a t = ^{1/{N — 1)) ~ (^)) 2 ? are deter- 

mined from the empirical data. The delayed exponential, f(t) = i](t — t )(l/r) exp(— (t — 
t )/r), where rj is the Heaviside step function, and t = minjshortest switching time, (t)—a t }, 
is used as the pdf. Once the waiting time to is determined, the maximum likelihood esti- 
mate for the lifetime is found as r = (t) — to- The slow- mode lifetime from this estimate 
is r = 13.1 ns, with a corresponding t = 7.8 ns. These results are shown in Fig. &[b) and 
Table [U with r for the slow mode close to other estimates that will be discussed below. In 
comparison, the empirical mean for the fast mode is (t) = 1.2 ns, with a t = 0.4 ns. The 
corresponding pdf is not well approximated by a delayed exponential for the fast mode. 

In addition to the above switches, 20 separate simulations were also completed that are 
identical to the previous conditions except for the details of the initial magnetization and 
field-reversal. The initial magnetization of these runs was determined by quenching the 
system to K while in the metastable well, and then re-thermalizing it to 20 K. From this 
time on (t = 0), the simulations were carried out identically to the previous 42 trials. The 
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reasons for this procedure are explained in Section IHI Al Along with the slow mode, the 
cumulative distribution of lifetimes for the quenched- relaxed (QR) runs is shown in Fig.[2](b), 
with r = 12.6 ns, a t = 10.8 ns, and t = 1.6 ns. 

A. Phase Portraits 

Phase plots of the energy also provide information about the behavior of the simulated 
nanopillar system. Such plots are shown in Fig. [31 with the energies due to dipolar and 
Zeeman contributions on separate axes. The collection of all runs belonging to a particular 
mode are shown in the background of each plot, with a single run overlaid on top. For 
the fast (a) and slow (b) modes, the simulations begin near the top of the plots, where the 
density of points is low. During the initial relaxation of these simulations, both modes evolve 
down and to the right (i.e., decrease in Zeeman energy and increase in dipolar energy). The 
obvious difference between plots (a) and (b) is the path that each takes near the metastable 
well (E z « —1130 erg/cm 3 ) and saddle point (E z « —1350 erg/cm 3 ), determined projective 
dynamics in Section [III CI The slow- mode trajectories proceed to the metastable well, which 
can be seen in plot (b) as the large dark region stretching from the center to the top left 
of the plot, but is absent in plot (a). For the slow mode, the simulation spends most of its 
time here. However, the fast-mode events ignore this attractor almost completely, mostly 
slowing down only near the saddle point. This is not unreasonable, since the free energy 
near the saddle point necessarily has a small gradient, and the driving force is therefore 
weak. Both modes continue the switching process toward the global free-energy minimum 
located below the displayed portion of the phase portraits. These results suggest that the 
difference between the fast and slow modes is the visitation of the metastable well by the 
slow mode. This is further supported by the QR simulations discussed below. 

The QR system begins at a point located more negative than E D rs —5200 erg/cm 3 along 
the dipolar axis of Fig. H^c) (not visible), which represents the common initial configuration 
for all QR trials. This initial configuration was found by quenching several simulations 
belonging to both the slow and fast modes at various values of Ez during the evolution of 
the system. Quenches initiated at values higher than the saddle point in Ez equilibriated 
to a common T = K metastable configuration. This configuration quickly proceeds to the 
T = 20 K metastable well when thermalized. Other quenches, which were initiated at values 
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FIG. 3: (Color online.) Phase plots in the space of dipolar and Zeeman energies for 29 runs 
belonging to the fast mode (a), thirteen runs belonging to the slow mode (b), and 20 runs belonging 
to the quenched-relaxed trials (c). The lighter background is the collection of all runs belonging 
to a particular mode, while the darker path represents a single run. Arrows indicate the average 
direction of motion of the phase portrait. 

of Ez below the saddle point, settled into configurations near the T = 20 K final absorbing 
state. 

Since this system necessarily starts in the metastable well, it does not have the same 
behavior as the fast and slow modes during the re-thermalization. However, ignoring this 
initial relaxation, it closely resembles the phase plot for the slow mode. The large change 
with T in the phase-space location of the metastable state of the QR systems is interesting 
in itself. It indicates that the entropy S, which enters the free energy F as F = E — TS, has 
a large influence on the free energy of the system. This observation is further supported by 
the fact that the system remains at a constant total energy E (except for small fluctuations) 
until the saddle point is crossed, as seen in Fig. HJ 

B. Transition Matrix 



Although the phase space in these simulations is very large (2N dimensions), it appears 
from the phase plots in Fig. [3] that the evolution of the system along the Zeeman-energy 
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FIG. 4: (Color online.) Total energy E as a function of time for a single run belonging to the slow 
mode and a single run belonging to the fast mode. Near t = 0, the total energy quickly decreases 
as the system initially relaxes. After this, the average value of the total energy remains constant, 
except for small fluctuations, until the saddle point is crossed. 

coordinate, E z , can approximately describe the process of magnetization switching. To 
investigate the possibility of a one-dimensional description of the switching process, this 
coordinate is discretized into 200 equal-sized bins labeled by i = 1, We found the 

results of the following sections to be approximately independent of the axis discretization, 
and we therefore use this value since it provides numerically stable results for all of the 
Markov-chain techniques. Each measurement of E z during the simulation corresponds to a 
particular bin i, and this discretized state can be represented by a unit vector (i\ consisting 
of 1 in the ith position with all other elements equal to 0. 

The matrix, M, of transition probabilities between states is also constructed by sampling 
the series of E 7, values during the simulation. Individual elements of the transition matrix, 
Mij, are obtained by enumerating the single transitions from bin j to bin i corresponding to 
a time step, At, equal to the measurement interval during the simulations, and normalizing 
such that Ylj Mij = 1. Thus, the probability of going from (z| to (j\ in one time step is 
given by the matrix element with the result that, 



(u(t + M)\ = (u(t)\M, 



(3) 



where (u(t) \ = = Si Pi ^ s ^ ne row vec t° r representing the probability that 
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E z is in bin % at time t. Therefore, M provides the average change in the state of the system 
after one time step.- These transition probabilities are found by combining statistics from 
all individual runs belonging to a particular mode (fast, slow, or QR). 

Since the average time evolution of the simulation decreases E z , the system begins in 
state (k\ and ends when the absorbing state (1| is reached and m z < (the z-component of 
the total magnetization is < 0). From this condition, the matrix element M\\ = 1 signifies 
this absorbing state. The transition matrix representing this absorbing Markov chain thus 
has the form, 4 



I Irxr r x,3 | 

M (r+s)x(r+s) =| | (4) 

Rsxr T sxs 

Here, I is an identity matrix which, in general, represents r absorbing states (here, r = 1). 
R is the recurrent matrix, which describes the probability of moving into the absorbing state 
from any other state, is a null matrix, and T is the transient matrix, which describes the 
evolution of the system before absorption. 

By this construction, M is a non-symmetric, regular, non-negative square matrix with 
different left and right eigenspaces. The Perron-Frobenius theorem 4 - provides a general 
property for this type of matrix, namely that there exists a unique eigenvalue of M equal 
to unity, which is larger than the magnitudes of all other eigenvalues of M. For our matrix, 
the left eigenvector associated with the largest eigenvalue of M is (Ai| = (1| = (1,0,0, ...0). 
This represents the probability distribution characterizing the absorbing state. 

The other eigenvalues of M correspond to decaying deviations from the equilibrium dis- 
tribution since (Ai|M = (Ax | and 

(X a \M = X a (X a \ (5) 

with < |A Q | < 1 implying that the average probability of being in each non-absorbing bin 
decreases with each simulation step. Here (Ai| = (1| denotes the dominant left eigenvector 
of M, and (A a | is any other eigenvector of M with eigenvalue X a . As a consequence of the 
above, all left eigenvectors of M except (Ai| will decay to zero under repeated applications 
of M since {X a \M. n = X a n (X a \. Equations (131) and ([5]) together are used to find the lifetime 
r of the ath eigenstate, with the result that, 

(X a (t)\ = Xj(X a (0)\ « e-^ x ^(X a (0)\. (6) 
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FIG. 5: (Color online.) Probability distributions obtained from the eigenvectors corresponding to 
the second largest eigenvalue for the fast and slow modes and the QR simulations. In this figure, 
the Zeeman axis has been discretized into 200 bins to be consistent with the projective dynamics 
analysis of Sec. IIII CI We found the probability densities to be approximately independent of the 
discretization. These distributions are constructed from all the data belonging to each mode, as 
discussed in Section fill Bl The slow and QR distributions have a well-defined peak at the location of 
the metastable free-energy well. The fast mode, however, is quite wide with its greatest probability 
near the saddle point. 

Since (A 2 | corresponds to the second longest-lived state of the system, r EV = 1/(1 — A 2 ) rep- 
resents the average lifetime of the metastable state of the system in units of the measurement 
time resolution. 

A probability distribution representing the metastable state can be built from a properly 
normalized linear combination of (Ai| and (A2I, (meta| = (Ai| + £»( A2 1 • Here, the scaling 
constant b constrains the resulting vector to have zero weight in the absorbing state and 
ensures that the probability is normalized (meta|z) = 1). As can be seen in Fig. [U the 
slow and QR metastable state probability distributions have a similar shape and exhibit 
peaks at the location of the metastable free-energy well. The fast mode, however, does not 
have a dominant peak in the probability distribution. Rather, the probability density is 
close to zero near the metastable free-energy well and spread out along the remainder of 
the axis, exhibiting a marginally larger density near the location of the free-energy saddle 
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TABLE I: Lifetimes in nanoseconds were obtained by the following methods: (t) is the empirical 
mean obtained from the switching times of the actual simulations, and at is the corresponding 
standard deviation. For the slow and QR modes, characterized by the function f(t) = rj(t — 
to)(l/r) exp(— (t — to)/r), to is found from minjshortest switching time, (t) — at}, and r = (t) — to- 
Since the fast mode is not well described by this function, to an d r are not defined for this mode. 
tev is obtained from the eigenvalues of the transition matrix, trt is found from the residence 
times of the PD analysis. 

Fast Mode (ns) Slow Mode (ns) Quenched Relaxed (ns) 



(t) 1.2 20.9 14.2 

a t 0.4 13.1 10.8 

to N/A 7.8 1.6 

r N/A 13.1 12.6 

r EV 0.5 21.6 13.5 

trt 1.2 22.7 15.2 



point. These results reinforce the conclusion that the fast mode switches simply do not fall 
into the metastable well. 

A summary of the lifetimes obtained by the methods discussed above is provided in 
Table [H We expect some bias in the lifetime data due to the artificial cut made when 
sorting fast and slow modes. This probably accounts for the large to for the slow mode. 
Trials with lifetimes < 2 ns are considered fast modes, although there is some probability 
that they may belong to the slow mode, albeit with a short individual lifetime. This may be 
reflected in the differences of the lifetimes of the slow mode and the QR trials. Differences 
between tev and the other measurements of the fast-mode lifetimes may be due to statistical 
error which is not included in the transition matrix associated with tevj or that for the fast 
mode analyzing only a single eigenmode may not be sufficient. We expect the slow and QR 
lifetimes obtained from the eigenvalues of the transition matrix to be more robust against 
this error due to the larger amount of statistics gathered from longer runs. 
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C. Projective Dynamics 

The transition matrix of the previous section accounts for all transitions that may oc- 
cur from a given bin, potentially with all elements of the submatrices R and T non-zero. 
However, since the bin sizes are chosen sufficiently large, the transition matrix is tridiago- 
nal. Under this condition, the one-dimensional description of the magnetization switching 
becomes a one-step Markov process.- Using this single, coarse-grained variable, the projec- 
tive dynamics (PD) method 3 - 1 ^ 1 ^ can be used to measure the growth probabilities Pq and 
shrinkage probabilities Ps of the stable phase along this coordinate. Once obtained, several 
properties of the projected free energy F(E Z ) can be measured via Pq an d Ps- 

The PD method is implemented as follows. First, the E z axis is broken into a number 
of bins, as described for the transition-matrix approach in Sec. IIII Bl The bin size is deter- 
mined such that each time step is only capable of moving the system between adjacent bins. 
Histograms of Pq(E z ), Ps(P z ), and the probability to stay in the same bin P^(E Z ), keep 
record of changes along the axis and are updated at each measurement. Once completed, 
the histogram is then normalized so that Pq(E z ) + Ps(P z ) + Pn(P z ) = 1. 

For the present nanopillar simulations, we attempted to use as the slow (binned) variables 
the Cartesian components of the total magnetization, as well as other contributions to the 
total energy. However, the Zeeman energy provided the best results for the projective- 
dynamics technique. This is reasonable since it provides the closest correspondence to the 
path observed in Fig. El 

Points along E z where Pq = Ps define local extrema of the free energy. This is true since 
P G > P s implies dF/dE z > and P s > P G implies dF/dE z < 0. For our system, which has 
a single metastable free-energy well, these extrema represent the location of the metastable 
well, the saddle point, and the true equilibrium in the free energy. The latter is not observed 
in our simulations due to the cut-off at m z < 0. 

Figures [6] (a) and (b) show PD plots for the fast and slow modes, respectively. Each plot 
contains the data of all the runs belonging to each mode (points), which have then been 
smoothed using a five-point running average (solid curves). The location of the metastable 
free-energy well for the slow-mode switches (first crossing from the right in Fig. [6] (b), E z m 
— 1130 erg/cm 3 ) coincides with the peaks present in the metastable probability distribution, 
Fig. El of Sec. IIII Bl Further left, the second crossing in Fig. [6] (b) indicates the saddle 
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FIG. 6: (Color online.) Projective Dynamics results for the two modes, fast (a) and slow (b). The 
slow mode exhibits clear crossings of the growth and shrinkage probabilities, indicating extrema in 
the free energy, corresponding to the metastable well (right) and the saddle point (left). Conversely, 
the probabilities for the fast mode only appear to overlap in this same region. The solid lines are 
five-point running averages. 

point in the free energy, which is located at E z « —1350 erg/cm 3 The locations of both the 
metastable well (first crossing from the right) and the saddle point (second crossing from 
the right) are obvious for the slow mode. This is expected since the results of Sec. IIIIBI 
indicate true metastable behavior. 

From our data, the fast mode has an overlap of Pq and P$ in the range from E z = —1350 
erg/cm 3 to E z = —1250 erg/cm 3 . It is not possible to determine if crossings occur in this 
region. This could be because crossings exist but are obscured by the small amount of 
statistics available, since less time is spent in this region than for the slow mode. It could 
also be that no true crossings occur. In this case fast swithches proceed through a region 
of the free-energy landscape with a small, but negative slope. Along with the results of the 
metastable probability distributions, this indicates the absence of a metastable state in the 
fast mode. These results agree with the probability distributions (Fig. [5]) obtained from the 
transition matrices of the previous section. 

The PD results of the quenched-relaxed system can be compared to the slow modes. As 
seen in Fig. [7J the QR system and the slow-mode probabilities near the metastable minimum 




15 




— Growth 

— Shrinkage 



000 -1800 -1600 -1400 -1200 

3 

Zeeman energy [ergs/cm ] 



-1000 



FIG. 7: (Color online.) Projective Dynamics results for the two modes, slow and the quenched- 
relaxed (QR). It is easy to see that these simulations share the same statistical behavior in this 
region, including nearly identical crossing values (which are locations of the extrema of the free 
energy). As with Fig. [6], the solid lines represent five-point running averages. 

cross at nearly identical values of Ez- Both PD and the metastable probability distributions 
imply that, once in the metastable well (the slow and QR decays), the system must be 
thermally activated to overcome the free-energy saddle point. However, if the metastable 
well is avoided (the fast mode), the system only has to traverse a relatively flat free-energy 
landscape to switch. 

Lifetimes can be calculated easily from the results of the PD technique.- The residence 
time for each bin h(i) is defined as the average time spent in state i. In addition, Pq{i)\i{i) 
is the average number of times the system moves from state % to i — 1 (the absorbing state 
is state 1). If the system is to reach the absorbing state exactly once then 



P G {i)h(i)-P s {i-l)h{i-l) = l 



(7) 



must be true. Since the absorbing state has no shrinking probability, Ps(l) = 0, the residence 
time of the previous state can be found by h{2) = 1/Pq(2). This leaves an iterative solution 
for the remaining states, 

l + P s (»-l)ft(»-l) 

Kv = ftt?\ • ( 8 ) 



Consequently, the average lifetime of the process is the sum of the residence times, trt 
h(i). The lifetimes obtained from this method are included in Table HI 



16 



IV. CONCLUSIONS 



We have studied the magnetization-switching properties of a simulated iron nanopillar 
motivated by experimental research.-^ Under the realistic physical conditions described 
here, we found the presence of more than one characteristic switching time, which were 
labeled as "fast" and "slow" modes. Through phase portraits and numerical results provided 
by transition matrices and the projective dynamics method, differences between these two 
modes were identified. Our results indicate that the fast mode is associated with switching 
dynamics that do not carry the system through a metastable well. This idea is supported 
by data from separate simulations in which the system was quenched to K while near 
the metastable well and then rethermalized. The quenched-relaxed simulations also indicate 
that the entropy provides a large contribution to the free energy of the system. 

Using transition matrices obtained from the average behavior of each mode, we also 
constructed projective dynamics plots and metastable probability distributions which further 
provide evidence that the fast mode does not have a metastable free-energy well. In addition, 
lifetimes were obtained by measurement, fitting to the cumulative distribution, sub-dominant 
eigenvalues, and residence times from projective dynamics and reported in Table [B These 
results support further evidence of the existence of multiple switching paths which may be 
verified experimentally, and may have important ramifications for technological applications 
that rely on a single, consistent switching mode such as memory devices. We hope to 
extend the scope of these results in future work by measuring the dependence of the fast- 
to slow-mode ratio as a function of the applied field angle and temperature. 
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